library(MCMCpack)

set.seed(200192)

filestring <- "uniform.prior.gamma."
M <- 50 ## number of Monte Carlo replications
## share of item pairs rated by a person among all the unique item pairs
share<-0.03
I.vec <- c(40, 80)       ## number of voters
J.vec <- c(40, 80)  ## number of votes
## If the corr bewteen the estimated par and the truth is less than this threshold,we store the data and the model output
corr.threshold<- 0.5 

# storage vectors
I.storage <- NULL
J.storage <- NULL
m.storage <- NULL
gamma.corr.storage <- NULL
theta.corr.storage <- NULL
gamma.mse.storage <- NULL
theta.mse.storage <- NULL
time.storage <- NULL

## constants across all Monte Carlo runs
burnin <- 500
mcmc   <- 2000
thin   <- 2
MCMCpaircompare2d_seed <- 3840
verbose<-mcmc+burnin+1 #shut the verbose
tune<-0.5 
procrustes<-TRUE




for (I in I.vec){
  for (J in J.vec){
    for (m in 1:M){
      start.time <- Sys.time()
      cat(filestring, "  I =", I, "  J =", J, "  m =", m, "\n")
      # generate true thetas
      theta.1.true <- rnorm(J)
      theta.2.true <- rnorm(J)
      # true values for anchoring items
      theta.1.true[1] <- 3
      theta.2.true[1] <- 0
      theta.1.true[2] <- 0
      theta.2.true[2] <- 2
      theta.1.true[3] <- -1
      theta.2.true[3] <- -2.5
      theta.true <- cbind(theta.1.true, theta.2.true)
      # generate true gammas
      gamma.true <- sort(runif(I)*pi/2)
      z.true <- cbind(cos(gamma.true), sin(gamma.true))
      
      ##create item pair comparison sparse matrix format
      npairs.i<-round(J*(J-1)/2*share) ## The number of pairs a rater judges
      rater <- NULL
      o1 <- NULL
      o2 <- NULL
      choice <- NULL
      for (i in 1:I){
        rater.id <- paste("rater", 10000 + i, sep="")
        for (p in 1:npairs.i){
          nums <- sample(1:J, size=2, replace=FALSE)
          o1.num <- nums[1]
          o2.num <- nums[2]
          eta <- t(z.true[i, ]) %*% (theta.true[o1.num, ] - theta.true[o2.num, ])
          prob.o1.chosen <- pnorm(eta)
          u <- runif(1)
          if (u <= prob.o1.chosen){
            chose.num <- o1.num
          }
          else{
            chose.num <- o2.num
          }
          rater <- c(rater, rater.id)
          o1 <- c(o1, paste("o", 10000 + o1.num, sep=""))
          o2 <- c(o2, paste("o", 10000 + o2.num, sep=""))
          choice <- c(choice, paste("o", 10000 + chose.num, sep=""))
        }
      }
      
      mydata <- data.frame(rater, o1, o2, choice)
      
      
      ## generate starting values 
      theta.start <- matrix(rnorm(2*J), nrow = J, ncol = 2)
      gamma.start <- runif(I)*pi/2
     
      
      out <- MCMCpaircompare2d(mydata,
                               theta.constraints=list(o10001=list(1,3),
                                                      o10001=list(2,0),
                                                      o10002=list(1,0),
                                                      o10002=list(2,2),
                                                      o10003=list(1, -1),
                                                      o10003=list(2,-2.5)
                               ),
                               verbose=verbose, burnin=burnin, mcmc=mcmc, thin=thin,
                               seed = MCMCpaircompare2d_seed, gamma.start=gamma.start,
                               theta.start=theta.start,
                               store.theta=TRUE, store.gamma=TRUE,
                               tune=tune, procrustes=procrustes) 
      # poterior sample means for gammas and thetas
      gamma <- colMeans(out[, grep("gamma", colnames(out))])
      theta <- colMeans(out[, grep("theta", colnames(out))])

      # correlation between the estimated parameters and the truth
      gamma.corr<-cor(gamma, gamma.true)
      theta.corr<-cor(theta, as.vector(theta.true))
      
      # mean square errors of the estimated parameters
      gamma.mse<-mean((gamma- gamma.true)^2)
      theta.mse<-mean((theta- as.vector(theta.true))^2)
      
      ## store info
      I.storage <- c(I.storage, I)
      J.storage <- c(J.storage, J)
      m.storage <- c(m.storage, m)
      gamma.corr.storage <- c(gamma.corr.storage, gamma.corr)
      theta.corr.storage <- c(theta.corr.storage, theta.corr)
      gamma.mse.storage <- c(gamma.mse.storage,gamma.mse)
      theta.mse.storage <- c(theta.mse.storage, theta.mse)
      end.time <- Sys.time()
      tot.time <- end.time - start.time
      time.storage <- c(time.storage, tot.time)
      print(tot.time)
      
      
      if(gamma.corr < corr.threshold | theta.corr<corr.threshold){
        filename <- paste(filestring, "unif.true.theta.I", I, "J", J,
                          "m", m, "Rda", sep=".")
        save(theta.true, file=filename)
        filename <- paste(filestring, "unif.true.gamma.I", I, "J", J,
                          "m", m, "Rda", sep=".")
        save(gamma.true, file=filename)
        filename <- paste(filestring, "unif.out.I", I, "J", J,
                          "m", m, "Rda", sep=".")
        save(out, file=filename)
        filename <- paste(filestring, "unif.dat.I", I, "J", J,
                          "m", m, "Rda", sep=".")
        save(mydata, file=filename)     
      }
      
                 
      
      
    } ## end M loop
  } ## end J loop
} ## end I loop


uniform.prior.gamma.out <- data.frame(I=I.storage,
                     J=J.storage,
                     m=m.storage,
                     gamma.corr=gamma.corr.storage,
                     theta.corr=theta.corr.storage,
                     gamma.mse=gamma.mse.storage,
                     theta.mse=theta.mse.storage, 
                     time=time.storage)


save(uniform.prior.gamma.out, file=paste(filestring, "out.Rda", sep=""))



